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Abstract. Recent work on synthetic rescues has shown that the targeted deletion 
of specific metabolic genes can often be used to rescue otherwise non-viable mutants. 
This raises a fundamental biophysical question: to what extent can the whole-cell 
behavior of a large metabolic network be controlled by constraining the flux of one or 
few reactions in the network? This touches upon the issue of the number of degrees of 
freedom comprised by one such network. Using the metabolic network of Escherichia 
coli as a model system, here we address this question theoretically by exploring not 
only reaction deletions but also a continuous of all possible reaction expression levels. 
We show that the behavior of the metabolic network can be largely manipulated by the 
pinned expression of a single reaction. In particular, a relevant fraction of the metabolic 
reactions exhibit canalizing interactions, in that the specification of one reaction flux 
determines cellular growth as well as the fluxes of most other reactions in optimal 
steady states. The activity of individual reaction scan thus be used as surrogates 
to monitor and possibly control cellular growth and other whole-cell behaviors. In 
addition to its implications for the study of control processes, our methodology provides 
a new window to study how the integrated dynamics of the entire metabolic network 
emerges from the coordinated behavior of its component parts. 
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1. Introduction 

Complex systems are composed of a large number of interacting parts. Physically, 
this means that the holistic description of a complex system necessarily involves a high- 
dimensional phase space. Significant previous work on the structure of complex networks 
[U El [31 m E] has been based on this paradigm. Herewith we explore an alternative 
approach, which we argue is appropriate to address the functional behavior of complex 
biological systems. By focusing on cellular metabolic networks as a model system of 
broad significance [6l [71 [H [9l [TO], [Til [12] , we show that despite being high dimensional, 
a metabolic network has a surprisingly small effective number of degrees of freedom, 
operationally defined as the number of independent reaction fluxes under steady-state 
conditions. This indicates that the steady-state dynamics of large complex networks 
can be significantly more constrained than their structure may suggest. Because the 
interactions constraining the dynamics also mediate information flow across the network, 
this result raises the possibility of natural as well as engineered control mechanisms based 
on the monitoring or manipulating one or few metabolic reactions. 

This study is motivated by the recent discovery that in single-cell organisms 
the growth defect caused by the deletion of an enzyme-coding gene can often be 
compensated by the concurrent deletion of a second enzyme-coding gene [12]. Such 
synthetic rescue interactions, in which damage compensates for damage, were predicted 
to even turn some non-viable gene-deficient strains into viable strains after specific 
additional gene deletions were introduced. Related research has found that cells evolved 
to optimize growth rate or any other typical function of metabolic fluxes tend to 
significantly reduce the number of active metabolic reactions when compared to typical 
non-optimal cells [T3|. This spontaneous reaction inactivation explains the role of latent 
pathway activation and why, sometimes, "less is more" in cellular metabolism. In 
particular, it shows that the compensatory perturbations underlying synthetic rescues 
are generally generated by the inactivation of metabolic reactions that are predicted 
to be inactive in growth-maximizing states. An outstanding question that stems from 
these findings is the extent to which the predetermined "optimal state" activity of a 
small fraction of reactions can constrain the entire metabolic network to operate close 
to the corresponding optimal state. Here we investigate this question by considering the 
flux specification of a single reaction under steady-state conditions. 

We focus on the latest reconstructed metabolic network of the bacterium 
Escherichia coli K12 MG1655 [H], which is arguably the most complete in silico 
cellular network to date [15], and we use flux balance-based methods [Ml [17] to 
computationally predict functional states of the system (see Appendix A). Within this 
framework, we observe that the pinned flux specification of a single reaction, such as 
the aminodeoxychorismate lyase reaction, can be sufficient to confine the steady-state 
cellular growth to zero or the maximum possible, without entailing other assumptions. 
In particular, we predict that the deletion effect of various otherwise essential enzymes, 
such as enolase (ENO), which catalyzes the conversion of 2-phosphoglycerate into 
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biomass cluster 




(CYSTL) Slave Reaction 

Figure 1. Master-slave relations, (a) Example of type A relation: fixing flux vi 
(master) at any feasible value uniquely determines flux V2 (slave), and vice versa, (b) 
Example of type B relation: fixing flux vi (master) at the extreme value vl uniquely 
determines flux V2 = V2 (slave), but the converse does not hold true. The reactions are 
indicated in parenthesis using the iAF1260 database abbreviation [M]. (c) Clustering 
characteristics of related fluxes. The diagonal blocks indicate reactions in type A 
relations, having either nonzero (blue) or zero (red) fluxes in any growth-maximizing 
state. The off-diagonal blocks (orange) indicate reactions in type B relations, whose 
fluxes are determined by fixing any of the corresponding block-diagonal reactions at 
the extreme values of a growth-maximizing state given by the simplex algorithm. The 
topmost blue block is the biomass cluster, which consists of the biomass reaction 
itself and 73 other type A-related reactions that uniquely determine the growth rate. 
Diagonal blocks with less than 4 reactions, which include all fluxes that can be both 
zero and nonzero at growth-maximizing states, are not shown. Additional information 
is provided in Table SI. 

phosphoenolpyruvate, can be compensated (up to a theoretical limit) by the controlled 
over-expression of a different reaction in the networl|j|. To explain and expand on these 
observations, we first determine how the steady-state condition alone affords such a 
canalizing interaction between the flux of one reaction and the integrated biomass flux 
by modeling master-slave relations between metabolic reactions. 

2. Flux relations 

We identify two fundamental types of master-slave relations (see figure [1]). We refer 
to them as type A relation, in which fluxes Vi and Vj can be uniquely determined by 
fixing either one of them at any feasible value; and type B relation, in which flux Vj 

X This often requires the coordinated activity of other reactions, as explained in the text. 
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is uniquely determined in an optimal state where flux Vi is maximized or minimized. 
These relations define clusters of fully and conditionally coupled fiuxes under steady- 
state conditions, as shown in figure [D^c). This figure shows that a total of 73 reactions 
are in type A relation with the biomass fiux, forming a cluster that we refer to as the 
biomass cluster. Notably, 92% of the other reactions are in type B relation with this 
cluster for any growth-maximizing state, i.e., they are "slaved" by the cluster. This 
is significant given that the set of growth-maximizing states itself is high dimensional 
[T8t fT9] . We also identify a number of smaller clusters. The second largest cluster 
consists of 22 reactions involved in lipopolysaccharide biosynthesis/recycling and cell 
envelope biosynthesis, which are nevertheless inactive in all growth-maximizing states. 
These results are ultimately related to the previously discussed concepts of fiux couplings 
[20l [2T] and fiux correlations [HI [22l [23l El], and thus potentially relevant for the 
identification of alternatives to known drug targets [25] and for the prediction of whole- 
cell metabolic behaviors based on the activity of a small set of reactions [26]. But what 
are the underlying mechanisms and functional consequences of this striking canalizing 
structure? 

Physically, the canalizing interactions can be interpreted as a consequence of the 
steady-state condition: the specification of a certain reaction fiux sets constraints on 
the possible fiuxes of other reactions in order to prevent accumulation and depletion of 
its products and reactants. At the most fundamental level, this is determined by the 
physical capability of the network to provide alternative pathways to produce and/or 
consume these compounds. In terms of the region of feasible fiux solutions, which is a 
convex region in the space of fiuxes determined by f™"^ and v^^^ {Appendix A), type A 
relations are satisfied in the interior of this region and type B relations are satisfied at 
the borders (cf figures 1(a) and (b)). 

Figure [2] shows a network representation of the biomass cluster, which is distributed 
across five functional subsystems of metabolism. Two intriguing properties emerge from 
this figure. First, despite being fully correlated, the reactions (blue squares) do not form 
a single connected network, but are instead separated into 16 different components. This 
property, which was also suggested in previous work [211 125] , can be readily rationalized 
when different components are linked together by parallel pathways whose combined 
fluxes are flxed. As shown in the simplifled representation of flgure[2]^b) for two reactions 
in the cofactor and prosthetic group biosynthesis, the individual reactions in these 
pathways need not to be locked to the activity of the biomass cluster. However, their 
combined activity is locked in steady states through the assumed constant concentration 
of metabolites Mi and M2. 

Second, while many reactions in the biomass cluster are coupled through the mass 
balance of metabolic compounds not involved in other reactions (green circles), we 
flnd a number of reactions coupled through metabolites that are also produced and/or 
consumed by reactions outside the cluster (red and orange circles). This counter- 
intuitive effect occurs when the reactions in the cluster and the other reactions sharing 
a common metabolite satisfy mass-balance conditions independent of each other. One 
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Figure 2. Biomass cluster, (a) Blue boxes and color circles represent the reactions 
in the biomass cluster and the associated metabolic compounds, respectively. Green 
circles indicate the metabolites exclusively produced and consumed by reactions in 
the cluster, while orange and red circles indicate that the metabolite is also produced 
and/or consumed by other reactions. Reactions and metabolites shown in gray do 
not belong to the cluster, (b-d) Examples of: (b) reactions in the cluster, such as 
AMPMS2 and TMPPP, that are coupled through multiple parallel pathways (gray 
symbols); (c) a metabolite (orange circle), such as fmn, that is shared by a decoupled 
reaction loop having complete mass balance of the metabolite within the loop; and 
(d) a reaction (yellow box), such as PDX5P02, that connects metabolites (red circles) 
directly linked to the biomass cluster and that is fully coupled to the cluster only 
in growth-maximizing states. The biomass cluster can be augmented to include 4 
additional such reactions (yellow boxes in (a)). Additional information is provided in 
Table S2. 
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such case is shown in figure [2]^c) for the local network of a shared metabolite Mi, such 
as fiavin mononucleotide (fmn), where the outside reactions are constrained to produce 
and consume Mi at exactly the same rate in order to balance metabolite M2. This 
mechanism can be conceptually understood as the decoupling between two different 
elementary flux modes [271 EH] • Figure E^d) illustrates a different structure in which a 
reaction Ri, such as pyridoxine 5'-phosphate oxidase (PDX5P02), is (bidirectionally) 
coupled to the biomass only when biomass production is maximized. In this case the 
flux of the parallel reaction R2, which is irreversible, goes to zero in growth-maximizing 
states, although it is generally nonzero in other states. 

This should be compared with parallel irreversible reactions connecting end-points 
of the biomass cluster, such as DMPPS and IPDPS in flgure[2](a). While the fluxes of the 
parallel reactions are not flxed by the biomass production rate in any steady state, their 
combined flux is, and thus the maximization of the individual fluxes of either DMPPS 
or IPDPS constrains all the fluxes of the biomass cluster. In the terminology developed 
above, the cluster is in type B relation with these reactions, i.e., the coupling direction 
is the opposite of the one shown in flgure[T](c). 

3. Pinned reaction expression 

We now turn to the implications of these mechanisms for the enhancement of biomass 
production (or growth) in reaction-deflcient mutants of E. coli, which are identical to 
the wild-type (WT) cells except that one reaction flux is constrained to zero. While the 
deletion of any reaction in the biomass cluster is lethal, since it forces the biomass flux 
to zerc[§|, controlling the flux of one such reaction can constrain the system to optimal 
states. This is achieved without requiring the speciflcation of the cell's response to 
perturbations other than the assumption that the reaction flux can be controlled and 
the post-perturbed flux distribution reaches a steady state. To extend this analysis to 
other reactions, we model the changes in cellular growth that may follow a reaction 
deletion and explore the pinned expression of a different reaction as a means to restore 
the growth or to identify surrogates for growth optimization in the reaction-deflcient 
strain (see flgure [3]) . 

Specifically, we test the hypothesis that growth can be restored if the fiux of 
the pinned reaction can be constrained to the optimal fiux value identified using fiux 
balance analysis (FBA) in a growth-maximizing state of the reaction-deficient mutant 
{Appendix A). In our simulations we assume that the pre-mutation organisms are in 
growth- maximizing states, as observed in adaptive evolution experiments [30]. We also 
assume that the early post-mutation state is governed by the minimization of metabolic 
adjustment (MOMA) hypothesis ^31], which has been shown to describe the fiux pattern 
of deletion mutants and assumes that the perturbed system goes to the closest available 

§ The biomass cluster is a subset of the previously identified essential reaction core [21] or lethality 
core [Uzgil B, unpublished], and it is related to the existence of a medium-independent activity reaction 
set ESI. 



Slave nodes and the controllability of metabolic networks 



7 




WT DEL DEL + PIN WT DEL DEL + PIN 




Pinned Reaction Pinned Reaction 



Figure 3. Pinned over- or under-expression of metabolic reactions in reaction- 
deficient mutants. Two implementations are illustrated, where the WT is assumed to 
be in a growth-maximizing state and the pinned reaction expression (PIN) is introduced 
(a) before or (b) after the network response to the reaction deletion (DEL). Aij 
indicates the growth rate change due to the pinned expression of a reaction j following 
the deletion of a reaction i, and Af '^'^ is the theoretical maximum of . (c-d) Color- 
coded normalized growth rate changes A = Aij / Af^^^ for the scenarios (a) and (b), 
respectively. On the horizontal and vertical axes we specify reactions whose pinning 
significantly increase growth or whose deletions are classified as lethal in experiments 
[32j [33l [Sll [35] , respectively. In (c), the data plotted corresponds to mutants having 
^max > of the optimal WT growth rate and to pinned expressions having A > 0.3 
for over-expressions (> 0.1 for under-expressions) for at least one mutant. In (d), the 
set of mutants and pinned reactions is the same as in (c) to allow direct comparison 
between the two implementation scenarios. A representation of (d) for data selected 
according to the same criteria used in (c) is shown in Fig. SI. Additional information 
is provided in Table S3. 

steady state in terms of Euclidean distances in the space of fluxes. MOMA-predicted 
fluxes generally correspond to suboptimal states, whose biomass production is lower 
than the theoretical maximum determined by FBA. 

We examine two possible scenarios, corresponding to two different experimental 
implementations. In the first scenario, we consider that the pinned reaction expression is 
introduced before the metabolic network responds to the reaction deletion (figure [3](a)). 
In the second scenario, we assume that the pinned reaction expression is implemented 
after the network has responded to the reaction deletion (figure [3](b)). The first 
case is modeled as a MOMA-predicted response to the combined perturbation of the 
reaction deletion and pinned expression (red arrow; figure [3l^ a)), while the second case 
is best modeled as a sequence of two MOMA-predicted responses (blue and red arrows; 
figure [3](b)). 
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Figures ^c) and (d) show the corresponding growth enhancement predicted for a 
selection of 38 reaction-deficient mutants. The biomass cluster is preserved for all single- 
reaction deletions in our simulations, indicating that under steady-state conditions 
the controlled expression of any of its reactions would tune the metabolic network 
of the deletion mutant to a growth- maximizing state. The same holds true for 11 
other reactions identified as directionally coupled to the biomass cluster for a growth- 
maximizing state identified by the simplex algorithm (i.e., the cluster is type B-related 
to them), namely EARlGOx, EAR160y, ASP03, ASP05, IPDDI, IPDPS, OPHHX, 
PDX5P02, MG2t3_2pp, TRPSl, and TRPS2 (see Table S3). Other reactions, such 
as DMPPS, can be directionally coupled to the biomass cluster for different growth 
maximizing state. Under the constraints imposed by constant biomass composition and 
steady-state reaction fluxes, the pinned expression of one such reaction confines the 
system to an optimal growth rate that depends neither on the MOMA modeling nor on 
the order of the perturbations (cf figures [3](c) and (d)). 

In addition, we identify reactions not related to the biomass cluster but whose flux 
pinning would significantly compensate the growth defect of the mutants. For example, 
while the suboptimal growth rate of the PDH-deficient mutant is predicted to be 0.28 
of the wild-type growth rate, pinning the flux of 3HAD121 increases the normalized 
growth rate to 0.65, which is 78% of the theoretical maximum (figure [3]^c)). The most 
significant growth recoveries are found to involve reaction over-expressions (43 reaction 
pairs with A > 0.3), but the positive impact of single-reaction under-expression [12] 
was also identified and corresponds to 10 reaction pairs with A > 0.1 in figure [3](c). 
Although the partial recoveries can vary with the implementation scenario, a statistically 
significant overlap is observed in figures [3](c) and (d): all full recoveries are the same 
and the growth changes are concurrently positive or concurrently negative for nearly 
78% of the other cases (see see Figs. S2). 

Figure [3] includes rescue counterparts for reaction deletions associated with genes 
identified as essential in deletion experiments [321 [SSj [SH [35]. For example, the over- 
expression of reactions of the cell envelope biosynthesis associated with fab genes is 
predicted to restore growth of several non-viable mutant^. 

We have focused on reaction activity in order to obtain general results entirely 
determined by the mass balance equations. However, as the above examples indicate, 
our approach also provides useful information about gene activity through known gene- 
enzyme-reaction relationships [36|. This is corroborated by the fact that 36% of the 
enzyme-coding genes in the iAF1260 reconstruction [1^ are in one-to-one relationship 
with metabolic reactions. Moreover, in many cases similar growth rescues can be 
obtained by pinning the expression of any one out of a large number of different reactions 
(figure [3]). From the perspective of metabolic control and the identification of alternative 
metabolic optimization targets, the latter indicates that one can choose to focus only on 
a subset of reactions with desirable properties (e.g., one-to-one reaction-gene relation or 

II Although our simulations do not predict zero suboptimal growth for all these mutants, the predicted 
growth rates tend to be significantly smaller than the wild-type growth rate. 
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availability of promoters and markers). 
4. Discussion 

The observation that the metabolic network may be controlled by a small number of 
degrees of freedom seems to be in accordance with the presence of global regulators 
[37] as well as the recent experimental observation that few (sometimes far-reaching) 
mutations can significantly increase the growth rate of E. coli strains subjected to 
adaptive evolution p8j. In interpreting our results in the context of metabolic 
engineering, one should of course not underestimate the experimental difficulties 
involved in the steady-state control of a reaction flux P, |39]. In particular, the 
reaction expression is not always correlated with gene expression and the availability of 
enzymes [40ll4T]. This is so partly because the balanced activity of a metabolic reaction 
may depend on the coordinated expression of multiple genes and may be influenced 
by post-transcriptional effects. While these issues fall outside the scope of this study, 
we observe however that signiflcant progress has been made in developing expression 
systems that can lead to tunable reaction expression patterns |^ B3]. These techniques 
can expand the applicability of our results as well as of the recently introduced OptReg 
platform [H], which is a versatile framework that exploits reaction down- and up- 
expression in the overproduction of targeted compounds. 

Speciflc experimental studies on E. coli metabolism appear to support our results. 
Examples of synthetic rescues induced by the suppression of speciflc metabolic reactions 
have been discussed in Ref. [I2], and similar experimental results are also found for 
reaction over-expressions. For instance, the growth of pgi-deficient E. coli mutants fed 
glucose has been shown to be signiflcantly improved by the over-expression of the soluble 
transhydrogenase UdhA |15]. This occurs presumably because UdhA restores redox 
balance when the Pentose Phosphate pathway becomes the primary route of glucose 
catabolism following the inactivation of phosphoglucose isomerase, which agrees with 
our predictions. Our modeling of the over-expression of the soluble transhydrogenase 
reaction leads to 5% flux increase through the Pentose Phosphate pathway while 
suppressing the flux through the Entner-Doudoroff pathway. 

Other experimental case studies that can be related to this work suggest that 
the possibility of rescuing a mutant using over-expressions is a general mechanism 
not limited to E. coli metabolism. For example, it has been shown that the over- 
expression of protein PGC-la, a regulator of energy metabolism, promotes the recovery 
of mitochondrial dysfunction caused by oxidant exposure ^6], apparently by up- 
regulating mitochondrial biogenesis in tissues with high metabolic demand. A different 
study has shown that HIV-1 mutants with budding defect are rescued by the over- 
expression of protein Nedd4-2s [47J, a Nedd4-like ubiquitin ligase of the family recruited 
by less-complex retroviruses. In humans, the over-expression of mitochondrial valyl 
tRNA synthetase has been shown to partially rescue cells carrying pathogenic mutations 
associated with inborn metabolic diseases [H]. Moreover, experiments with insulin 
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receptor-deficient mice, a mutant that develops severe diabetes, indicate that the over- 
expression of hepatic glucokinase improves glucose tolerance and partially compensates 
for the metabolic disorders associated with this deficiency [19]. In addition, a number 
of over-expression rescue interactions have been identified for genes involved in various 
other cellular functions, most noticeably for yeast [50], as shown in the Saccharomyces 
Genome Database ^51j • While the mechanisms underlying these examples remain largely 
unexplained, they highlight the potential significance of the interactions systematically 
predicted here, particularly for the recovery of lost cellular function and the substitution 
of known drug and microbial optimization targets. 

The reduced effective number of degrees of freedom identified in this study 
may represent a general property of complex biological networks whose function, 
like in metabolism, is based on the transport and/or transformation of locally 
preserved quantities. This includes, for example, food webs and many other resource 
transportation or transformation networks. In such systems, the suppression of one flux 
is often accompanied by the enhancement of different fluxes, allowing us to interpret 
flux over-expressions as reciprocals to the flux down-expressions recently exploited 
to bypass defective pathways [12j. We expect that the insights provided by this 
study will be further expanded in combination with other techniques, such as the 
extreme pathway [52] and elementary flux mode [271 [2H] analyses, through the study 
of non-stationary behavior [53], and by means of additional applications of network 
analysis [H [2], [3], HI [5], [551 ESI [571 ISHl [59]. Altogether, this promises to improve our 
understanding of the interaction between the dynamics and the control mechanisms 
underlying the functional behavior of complex biological networks. 

Acknowledgments 

The authors thank Natali Gulbahce and Linda Broadbelt for insightful discussions. This 
work was supported by NSF under Grant No. DMS-0709212 and involved the use of 
computer equipment funded through the NSF-MRSEC Program DMR-0520513. 

Appendix A. Network Model and Methods 

The genome-wide in silico E. coli metabolic network iAF1260 [H] used in this study 
consists of n = 2 074 unique biochemical reactions and m = 1 039 chemical compounds. 
The network and the state of the system are conveniently represented by a m x n 
matrix of stoichiometric coefficients S = (Sij) that indicates the molar proportions of 
the reactants i in reaction j and a n-dimensional vector of reaction fluxes v = (fj), 
respectively. Our steady-state analysis is based on identifying solutions for the reaction 
fluxes through the mass balance equation Sv = subjected to ff™ < < v^"^^, where 
f ™™ and t)™^^ represent limitations imposed by nutrient availability and thermodynamic, 
physiological or biochemical constraints [IE]. For concreteness, we consider a limited 
glucose nutrient environment with maximum uptake rates of 10 for glucose and 20 for 
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oxygen (in units of mmol/g DW-h). In glucose aerobic media, 152 regulated reactions 
of the reconstructed model are assumed to be inactive and a split of 1:1 is used for 
the flux ratio between the two NADH dehydrogenases [Hj. We focus on the set of 
= 1 287 biochemical reactions involving M = 690 chemical compounds that can be 
active under these conditions. Of this total, 853 reactions are necessarily inactive in 
growth-maximizing states. Specific metabolic states are determined through two widely 
used optimization schemes, flux balance analysis (FBA) for growth- maximizing states 
[T6] and the minimization of metabolic adjustment (MOMA) for post-perturbation 
states [31] . FBA finds a solution v that maximizes the production of biomass, which is 
represented through an additional reaction in matrix S, and exploits the fact that the 
biomass flux rate can be mapped to growth rate. MOMA finds a new feasible solution 
Vnew that is closest to the original state t^o in terms of the Euclidean distance | |i^ne«)— I'd! | 
in the space of fluxes. The master-slave relations are identified using a variant of flux 
variabihty analysis [18], which identifies the upper and lower bounds of the flux values 
when the flux of a given reaction is fixed, as in figure [Hj^a) and (b). This approach avoids 
the combinatorial explosion problem [51] inherent to extreme pathway analysis [52] and 
elementary flux mode analysis [271 EH] , allowing it to be applied to the full-scale network. 
The simulations are implemented using the ILOG CPLEX optimization software and 
the simplex algorithm. 



Appendix B. Supplementary Data | |60] 

Table SI. Reactions in the two largest diagonal blocks indicated in Fig. [T](c). 

Table S2. Additional information about the biomass cluster in Fig. [2|^a). 

Table S3. Deleted and pinned reactions in Figs. [3]^c) and[3t^d). 

Figure SI. Different representation of Fig. [^d). 

Figure S2. Direct test of the overlaps between Figs. [3](c) and[3](d). 

Test of different objective functions. ATP production and lactic acid production. 
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